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Abstract. We present iV-body simulations of unstable spiral 
modes in a dynamically cool collisionless disc. We show that 
spiral modes grow in a thin collisionless disk in accordance 
with the analytical perturbation theory. We use the particle- 
mesh code SUPERBOX with nested grids to follow the evolu- 
tion of unstable spirals that emerge from an unstable equilib- 
rium state. We use a large number of particles (up to = 
40 X 10^) and high-resolution spatial grids in our simulations 
(128"^ cells). These allow us to trace the dynamics of the unsta- 
ble spiral modes until their wave amplitudes are saturated due 
to nonlinear effects. In general, the results of our simulations 
are in agreement with the analytical predictions. The growth 
rate and the pattern speed of the most unstable bar-mode mea- 
sured in TV-body simulations agree with the linear analysis. 
However the parameters of secondary unstable modes are in 
lesser agreement because of the still limited resolution of our 
simulations. 

Key words: Stellar dynamics - Galaxies: kinematics and dy- 
namics - Galaxies: spiral - Galaxies: structure 
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J> 1. Introduction 
• I— I 

, Gravitational instability of galactic discs is a widely accepted 
physical mechanism for the generation of baiTed and spiral 
structures in disc galaxies. Since the pioneering publications by 
Lin & Shu (1964) and Toomre (1981), considerable progress 
has been made in the theoretical study of gravitational insta- 
biUties of disc-like systems. A steady stream of papers per- 
forming a linear stability analysis of self-gravitating gaseous 
discs (e.g., Bertin et al. 1989a,b; Adams et al. 1989; Savonije 
& Heemskerk 1990) have demonstrated that a massive gaseous 
disc will almost inevitably be prone to spiral instabilities. A 
number of numerical techniques have been applied to study the 
nonlinear evolution of unstable gaseous discs (e.g., Tomley et 
al. 1994; Miyama et al. 1994; Woodward et al. 1994; Nelson 
et al. 1998; Laughlin et al. 1997,1998). Recently, the linear ap- 
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proach has been generalized for the stability analysis of non- 
axisymmetric gaseous disks (Asghari & Jalali 2006). These 
contributions address the dynamical modeling of gaseous discs 
with a variety of equilibrium properties and compare the results 
of numerical simulations with those of analytical predictions. 

The dynamical study of the collisionless gravitating discs 
is hampered by the difficulty of normal mode calculations in 
such systems. A global modal analysis of collisionless discs has 
been accomplished for a few cases. Kalnajs (1972) calculated 
the global modes of a uniformly rotating, self-gravitating stellar 
disc. Zang (1976) carried out a modal analysis of Mestel's stel- 
lar disc using an inner cutout function to handle the central sin- 
gularity. His method was then used by Evans & Read (1998a,b) 
for a general class of scale-free discs with rising and falling 
rotation curves. The modal properties of a collisonless expo- 
nential disc (with a core) were also studied by Vauterin & De- 
jonghe (1996). Most recently, Jalali & Hunter (2005, hereafter 
JH) developed a rather general approach that allowed them to 
calculate the modal properties of stellar discs. They showed 
that the radial orbits of soft-centered models lead to a bound- 
ary integral that plays a crucial role in the formation and growth 
of a bar mode, which is localized within the central regions of 
a collisionless disc. They also investigated the effect of dark 
halos/bulges on the properties of growing modes. 

The major drawback of perturbation theories is that they 
cannot make any prediction of modal evolution in the nonlinear 
regime. A linear analysis cannot predict the saturation level of 
the wave amplitude after an exponential growth phase, nor can 
it predict the duration of the linear and nonlinear phases. One 
can tackle these problems by using iV-body simulations. 

An A^-body approach has a number of difficulties, such as 
numerical relaxation due to a high level of noise and finite 
spatial resolution. These factors hamper the behavior of col- 
lisionless models. Moreover, growing modes are sensitive to 
the equilibrium profiles of the central regions. Thus, numerical 
simulations with a large number of particles, and high spatial 
resolution, are needed to correctly model the behavior of un- 
stable modes in collisionless discs. 

There have been a few attempts to numerically model the 
dynamics of global modes in collisionless discs. As already 
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mentioned, modeling of the spiral structure in collisionless 
galaxies and comparison with analytical predictions has been 
restrained by the relatively small number of particles used 
in A^-body codes. Randomly generated initial coordinates of 
particles inevitably lead to density fluctuations and perturb- 
ing forces that affect the disc dynamics (Sellwood 1983). The 
disc response to these virtual forces could contribute substan- 
tial errors of about 50 percent to the computed growth rates 
and the pattern speeds (Sellwood 1983). In previous A^-body 
simulations aimed at studying the dynamics of global modes in 
collisionless discs (Athanassoula & Sellwood 1986; Sellwood 
& Athanassoula 1986), special precautions have been taken to 
avoid a strong influence of noise perturbations. These authors 
used 'quiet' initial conditions, with the groups of stars spaced 
at equal intervals around the rings at fixed radii, and all parti- 
cles in a ring were given the same initial radial and tangential 
velocity components. 

The development of numerical techniques with tens of mil- 
lions of particles would avoid artificial precautions in tackling 
the numerical noise problem. In this paper, we use SUPERBOX 
(FelUiauer et al. |2000| l to model the evolution of a collisionless 
disc. SUPERBOX is a highly efficient particle-mesh-code with 
nested and co-moving grids, based on a leap-frog scheme de- 
signed for the simulation of interacting galaxies or other stellar 
systems. Aiming at modeling the dynamics of spiral patterns in 
real galaxies, we apply the code to analyse one of JH models 
that has an exponential disc density distribution with a core, 
embedded in a cored logarithmic potential. We use this model 
for two reasons. Firstly, closed-form expressions are available 
for the phase space distribution functions. Secondly, a full spec- 
trum of unstable modes has been computed for this model. 

We find good agreement between the linear global modal 
analysis and the nonlinear simulations. The growth rate and the 
pattern speed of the most unstable m=2 bar-like global mode 
found in A^-body simulations, and the spatial distribution of 
this mode, agree well with the analytical results. We have also 
been able to recover the theoretical growth rates and the pattern 
speeds of the unstable m=3 and m=4 global modes. Apart from 
demonstrating the existence of global modes in the cored expo- 
nential discs, our simulations provide a welcome check of the 
SUPERBOX code, and demonstrate its applicability to model 
the dynamics of real galaxies. 

In fj2] we briefly summarize the properties of the global 
modes of JH models, and use their method to find the unsta- 
ble 771=2, m=3 and m=4 modes of a cored exponential disc. In 
0we describe the A^-body code and its results for the model 
of JH. A comparison of the results of A^-body simulations with 
those of analytical predictions is given in the same section. 21 
contains a summary of our results. 



2. The analytical model 

We study the stability properties of a razor-thin disc that has 
an exponential density distribution with a central core of radius 
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Fig. 1. Upper panel: Solid line shows the rotation curve Vmt 



of the model. Dotted line is the streaming velocity 



and 



dashed line is the normalized surface density distribution. 
Lower panel: Solid and dotted lines show Cr and G^, re- 
spectively. Dashed line is Toomre's Q minus 1 (Q — 1). We 
have set the model parameters to A^=6, CEqRd /vq=Q.34 and 
A = Rc/Rd=0.625. 
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We assume that the disc is in equilibrium with the total gravi- 
tational potential of the galaxy represented by the logarithmic 
law 



MR) = v'aln^l + RyRl. (2) 

Here $o {R) is the total potential of the disc and halo compo- 
nents. The disc rotation curve is given by the equation: 

VqR 



Vrot ^ \ R 



(3) 



dR ^Rl. + i?2 ■ 
Equation (O gives the rotational velocity of a collisionless disc 
if the velocity dispersion of the disc is equal to zero. In a disc 
with a non-zero velocity dispersion, the mean rotational veloc- 
ity (uy) differs from the simple law given by equation ([3]) due 
to a collisionless 'pressure' that influences the rotation of the 
disc. 

The density distribution given by equation ([T]i is charac- 
terized by three parameters, namely the central surface density 
So exp(— A), the radial scale length of the disc density distribu- 
tion Rd (through A), and the core radius Rc- We use the core 
radius Rc and the asymptotic circular velocity vo as units to 
normalize the problem. The time unit is then determined by the 
ratio Rc/vq. The dimensionless parameter 5*0 = GSo-Rd/wq 
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Fig. 2. Unstable modes of the cored exponential disc for A^=6, A = Rc / Rd=0-6^5 and CSqRd /vq=0.34-. Contours show the 
positive surface density perturbations. The contour levels are equally spaced from 10 to 90 percent of the maximum. Solid and 
dotted circles mark the CR and OLR circles respectively, (a) The fundamental m = 2 mode, (b) The secondary m = 2 (spiral) 
mode, (c) m — 3. (d) m = 4. 



gives the ratio of the mass of the disc to the total mass of the 
model within the radius R w 5Rc'- 



Following JH, we assume the distribution function of the 
disc particles in the form: 
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Here E and L are the total energy and the angular momentum 
of an individual star as 

1 . o 



E 



$o(i?), L^Rv^ 



(6) 



The family of distribution functions given by equations ^ and 
(|5]l depends on an integer parameter N, that controls the frac- 
tion of near-circular orbits, and thus controls the disc velocity 
dispersion. The radial and azimuthal velocity dispersions de- 
termined as 

Cb^ = = ^K) - K)' (7) 

are small for relatively high values of N . We select a mod- 
erately cold disc characterized by the parameters = 6, 
A — Rc / Rd=0-625 and 5o=0.34. The upper panel of Figure 
[T]shows the rotation curve Vi-ot determined by equation (|3]l, the 
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mean rotational velocity {v^p) and the normalized surface den- 
sity e^Y^niR) /^Q for this model. The lower panel of Figure 
[U shows Toomre's Q-parameter together with the radial and 
azimuthal velocity dispersions as a function of the radial dis- 
tance R. For the selected model, the Q-parameter is greater 
than unity everywhere in the disc, and the model is stable for 
m=0 perturbations. 

In the linear regime, the spiral density perturbation can be 
written in polar coordinates {R, (p) as: 

E™(i?,^,t) = i„^o(i?)e'("'^-"*\ (8) 
where 

u^mnp + is, i„^o(i?) = ^m(i?)e"'"(^). (9) 
Here Am{R), Om{R), and s are the amplitude, phase, pat- 
tern speed and growth rate of the mth global mode, respec- 
tively. 

Figure|2]displays the results of the linear global modal anal- 
ysis of our model. We find that the coUisionless disc is unsta- 
ble towards two, three and four-armed spirals. We have shown 
the contour plots of the unstable eigenmodes, together with the 
positions of the corotation (CR) and outer Lindblad (OLR) res- 
onances. Pattern speeds Hp and growth rates s of all unstable 
modes are collected in Table |2] 

Figure |2ji shows the fastest growing m=2 bar-like mode 
(2p in Table |2]l. The dimensionless pattern speed and growth 
rate for this mode are ftp = 0.768 and s = 0.642, respectively. 
Linear analysis also reveals another unstable m=2 global mode 
(FigurelJI?) that has a pattern speed of Hp = 0.443 and a growth 
rate of s = 0.119 (2s in Table|2]i. This secondary mode is more 
spatially extended than the fundamental bar-mode, and occu- 
pies the disc within a few core radii. The spatial distribution of 
amplitude of the bar-mode has a single maximum located ap- 
proximately at a distance of i? ~ 0.25 from the centre of the 
disc. The amplitude function of the secondary m=2 mode has 
three maxima of comparable amplitudes shifted approximately 
by 90 degrees with respect to each other. The secondary m=2 
mode has also been shown in Figure 9 of JH, and we reproduce 
it here for completeness. 

Figures ^ and ^ show the contour plots for the unsta- 
ble m=3 and m=4 modes found in the linear stability analysis. 
These modes occupy approximately the same region of the disc 
as the secondary m=2 mode does. These modes, however, grow 
faster than the secondary m=2 mode (see Table |2]i. Therefore, 
the dynamics of perturbations in the outer regions of the disc is 
governed mostly by the m=3 and m=4 perturbations. The pat- 
terns of TO = 3 and m — A instabilities freely extend up to the 
outer Lindblad radius. This shows that the corotation resonance 
has no influence on these modes. 



3. Numerical simulations 

3.1. Code and data analysis 

Our A^-body simulations are carried out with the help of the 
SUPERBOX code (Fellhauer et al. 120001 1. This is a highly ef- 
ficient particle-mesh-code with nested and comoving grids. 




-j I I ' I I I I I ' I I I I 
2 4 6 8 10 12 r, kpc 

Fig. 3. The rotational velocity and velocity dispersions of the 
disc as a function of the radial distance R. Solid lines are the 
disk equilibrium profiles calculated with help of the distribu- 
tion function (equation |4]i. The points correspond to the iV- 
body realizations of the initial equilibrium. Dotted line shows 
the rotation curve Vo{R) of equation Q. 



based on a leap-frog scheme with second order force calcula- 
tion. Nested grids, comoving with the center of mass, allow us 
to achieve high resolution in the central parts of the collision- 
less systems. The code has been successfully applied to study 
mergers of galaxies (Madejsky & Bien 1993 ), galaxy-satellite 
disruption (Klessen & Kroupa ll998l l and the orbital decay of 
satellite galaxies (Penarrubia et al. 120021 Just & Penarrubia 
2005 ). It was also used to study orbital evolution of a super- 
massive black hole in a galactic nucleus (Spinnato et al. l2003l l. 

We use SUPERBOX to simulate the dynamics of perturba- 
tions in a coUisionless disc, which is in rotational equilibrium 
under the influence of its self-gravity and the gravitational po- 
tential of an external rigid halo. The code is three dimensional. 
To study the dynamics of a two-dimensional gravitating disc, 
we confine the initial distribution of particles to the {x, y) -plane 
and set vertical velocities of the particles equal to zero. The dy- 
namics of the model was simulated using different sets of the 
numerical parameters (i.e., number of particles, grid resolution, 
etc.) listed in Table[T| The radial extent of the inner, intermedi- 
ate and the outer grid zones are determined by the parameters 
Ri, Rm and i?o, and the spatial resolution in each grid zone is 
determined for the inner grid zone by the parameter di. 

The SUPERBOX has a fixed time step, so to resolve the mo- 
tion of the particles in the central regions of the disc, we use a 
time step of dt = 0.05 Myr This is less than a typical crossing 
time of a grid cell in the inner regions of the disc (« 0.1 Myr). 

The disc is built with equal mass particles. To construct a 
2D coUisionless model in dynamical equilibrium, we use the 
distribution function given by equation (|4]), and set velocities 
of the particles in z-direction and their z-coordinates equal to 
zero {zj — Vz.j — 0). The gravitational potential of a rigid halo 
is assumed to be spherical and is calculated by subtracting the 
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Table 1. Parameters of the different models. All models have the same circular velocity = 220 km/ s and disc scale length 
Rj^ = 3 kpc. Model A with A Rc/Rd 1 and 5*0 = 0.133 is the stable reference model. Model B with A = 0.625 and 
a larger disc mass So = 0.329 corresponds to the unstable model in ^ dim and Rmax are the dimension and cutoff radius of 
the disc. Eq determines the surface density of the disc (see Eq.[rii and dt is the time step of the models. N gives the number of 
particles, the next column the number of grid cells in each grid. i?m and Ro are the radii of the inner, middle and outer grid. 
di gives the cell size of the inner grid cells. 



disc potential $ d from the total potential given by the equation 

To determine the macroscopic characteristics of the disc, as 
a function of radius, such as surface density T,{R), disc rotation 
(v^) and the radial and azimuthal velocity dispersions Cr and 
Cip, we divide the disc into rir = 100 rings of equal width and 
calculate the macroscopic values in each ring. Figure [3] shows 
the A^-body realization of the initial equilibrium profiles of the 
disc. In this Figure, the solid lines show the the rotation ve- 
locity and the velocity dispersions of the disc, calculated with 
equation (|4|i, and the points are the iV-body realization of disc 
equilibrium. As it can be judged from Figure[3] there is satisfac- 
tory agreement between the theoretical model and its Ai"-body 
realization except for the very central regions of the disc. 

To quantify the growth rates and the pattern speeds of the 
unstable modes, we calculate the Fourier components of the 
perturbed density for azimuthal wavenumbers m — 1, 2, • • • , 6. 
With point-like particles, the Fourier components are given by 
the real part of the expression 
Mi 



{Ri 



(10) 



where Si is the area of the ith ring element corresponding to 
Ri, Ni is the number of particles in the zth ring, and Mi is the 
total mass of particles that lie on Si. By defining the amplitude 

A„.(i?,) = Arn{R^) we get 

T,jn ^ Ajn{Ri) COS {mip - Om) ■ (H) 

Generally, the amplitude Am{R, t) and the phase 6m{R, t) of 
the modes depend on time and radius. By measuring these 
quantities, we can estimate the growth rate s and the pattern 



speed rip of the modes by fitting the perturbed quantities to the 
expressions 

A^{t) A™,oe'*, Ora{t) 6l„,o + mnpt. (12) 
If the growth rate and the pattern speed of a Fourier component 
are independent of radius, the perturbation is dominated by an 
eigenmode. 

To test the influence of the resolution effects and to deter- 
mine the long-term behaviour caused by numerical errors, we 
simulated the dynamics of a stable model with a low disc mass 
(model A in Table [T]i. For this model, the number of particles 
as well as the grid resolution has been varied. We find that af- 
ter adapting the particle distribution to a grid-based potential 
the disc stays in equilibrium. After 3 Gyr, the artificial heat- 
ing caused by the particle-mesh scattering is lower than 5% of 
the initial velocity dispersion which shows that SUPERB OX is 
intrinsically collisionless. 

To test the dynamics of a three-dimensional stellar disc, we 
expanded model A2 in the vertical direction to an isothermal 
slab with the scale height of h = lOOpc (model A3). The re- 
sult of this exercise is that the dynamical equilibrium is not 
affected by the thickening and that the stability of the disc is 
not destroyed. 



3.2. The unstable model 

In this section we discuss properties of the unstable modes 
found in numerical simulations, and compare the results to 
those obtained in a linear analysis. We find that the low- 
resolution model BO leads to a significant discrepancy between 
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Fig. 4. Temporal evolution of the surface density distribution for an unstable disc (model shown in Fig. 2). A rapidly growing 
rotating bar-mode emerges in the central regions of the disc at early stages of disc evolution. At later times, the more slowly 
growing three-armed spiral determines the dynamics of perturbations in the outer regions of the disc. The time t is given in Myr. 



N-body simulations and analytical predictions. We therefore 
discuss results obtained for the higher-resolution model Bl 
which is built with 13 million particles. Models B2-B8 differ 
from the model Bl in the numerical parameters, namely the 
outer cutoff radius, the time step, the number of particles, and 
the grid resolution, and serve to illustrate that parameter vari- 
ations do not affect the results of N-body simulations for our 
fiducial model. Model B9 follows the disc evolution three times 
longer than Bl with the same parameters to investigate the non- 
linear evolution of the unstable modes. Models B10-B12 are a 
series, where only the number of particles is changed for the de- 



termination of the N dependence of pattern speeds and growth 
rates. 

Figure |4] shows a time sequence for the contour plots of 
the perturbed disc surface density plotted up to the nonlinear 
saturation phase. Linear analysis of our model predicts the ex- 
istence of a few simultaneously growing unstable modes. The 
results of N-body simulations are in qualitative agreement with 
the linear analysis predictions. As expected, the early stages of 
disc evolution are governed by the fastest growing m=2 global 
bar-mode, developing in disc central regions. This mode is seen 
in Figure m starting from 25 Myr By 60 Myrs, the central bar- 
mode saturates, and the perturbations continue to grow in the 



A.V. Khoperskov et al.: Simulations of unstable modes 



outer disc regions. The fastest growing mode is confined to the 
central kiloparsec of the disc, and the dynamics of the outer 
disc regions is governed by more slowly growing spirals. The 
disc dynamics in outer regions is determined by a superposition 
of slowly growing two-armed, three-armed and four-armed spi- 
rals. These modes are shown in Figure |2l A^-body simulations 
confirm this finding. 

For quantitative comparison of growing global modes 
found in A^-body simulations with the analytically predicted re- 
sults, we calculate a Fourier decomposition of the density per- 
turbations. This allows us to determine the amplitude A„i{r, t) 
and phase Qmir, t) of an m-armed spiral mode as a function 
of time and of the disc radius, and to determine the growth rate 
and the pattern speed of each Fourier component as a func- 
tion of radius. Figure |5] shows the amplitude and the phase of 
Fourier modes m — 2, 3, 4 calculated at different radii in the 
disc. In the central regions, the dynamics of perturbations is 
governed by the fastest growing m=2 mode, which saturates at 
about 40 Myr. In the outer disc regions, the dynamics is dic- 
tated by the m=3 perturbation prevailing over its m=2, and 
TO=4 competitors. As can be judged from Figure |5] the spi- 
ral growth rates in the linear stage and their pattern speeds are 
fairly independent of radius, which indicates that perturbations 
are indeed the growing global modes. The growth rates s and 
pattern speeds fip for the primary mode m = 2 do not change 
value at the transition from the linear stage to the essential non- 
hnear stage (t » 1/s). 

This is additionally illustrated in Figure |6] showing the ra- 
dial dependence of the growth rates, pattern speeds (upper 
frames), and the amplitudes of m=2, 3, 4 Fourier components 
as a function of radius. The parameters have been calculated 
from A^-body simulations at time 40 Myr. The fastest grow- 
ing bar-perturbation grows as a whole in the central regions of 
the disc, and rotates with fairly constant angular velocity. The 
measured values for the growth rate and for the pattern speed 
for this mode agree within errors of measurement with the lin- 
ear analysis. The pattern speed of the bar-mode of flp = 0.78 
agrees with the theoretical value of ftp — 0.77, and the growth 
rates found in both linear analysis and A^-body simulations 
quantitatively agree as well. The amplitude of the bar-mode 
reaches its maximum at about 0.25 kpc (bottom frame of Fig- 
ure|6]l, and and then decreases with radius. Such a radial profile 
agrees qualitatively with the linear analysis predictions (solid 
line). A discrepancy between the linear analysis and numerical 
simulations in the outer regions of the disc is explained by the 
presence of a secondary m=2 global mode. 

The agreement between growth rates and pattern speeds for 
m=3, and m=4 modes measured in N-body simulations and 
those determined in linear analyses are less satisfactory. While 
pattern speeds agree reasonably well, the growth rates mea- 
sured in N-body simulations show a considerable mismatch 
with the values found in linear analyses. Being lower than the 
growth rate of the fastest growing bar-mode, the growth rates in 
N-body simulations are higher than the analytical predictions. 
In part, the discrepancy can be related to the still low resolution 
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mode 




analytical 




simulations 


m 


i Lp 


s 


Ren 


^OLR 


(ilp) (s) 


2p 


0.768 


0.642 


0.834 


2.077 


0.78 0.63 


2s 


0.443 


0.119 


2.042 


3.775 


0.61 0.41 


3 


0.597 


0.153 


1.345 


2.320 


0.62 0.33 


4 


0.653 


0.164 


1.159 


1.880 


0.57 0.24 



Table 2. Growth rates and pattern speeds of the different modes 
for model B12. The time unit is Rc /vq = 8.52 Myr. 



of our simulations which is not high enough to properly model 
multi-modal behaviour of the unstable disc. 

Figure |7] shows the dependence of the growth rate and the 
pattern speed of the principal bar-mode as a function of the 
number of particles used in N-body simulations. As can be seen 
from Figure [T] the parameters of the unstable mode depend 
strongly on resolution, i.e., number of particles used in N-body 
simulations. The experiment with 3 million particles gives a 
large discrepancy between the predicted parameters of the un- 
stable mode (solid lines in Figure IT) and the values measured 
in N-body simulations. A simulation with 13 million particles 
gives better agreement with theory, and one with 30 million 
particles gives a satisfactory agreement. Additionally, SUPER- 
BOX achieves best resolution in the disc's central regions, so 
the parameters of the centrally confined fastest growing bar- 
mode are in better agreement with the analytical predictions 
compared to the more slowly growing modes developing in the 
outer regions of the disc. An increase of particle resolution in 
N-body simulations should lead to a better agreement of the 
parameters of the secondary unstable modes with theory. 

The spatial appearance of the eigenmodes is determined by 
the radial dependence of the amplitude and the phase of each 
Fourier component. 

The 771=2 Fourier component is a superposition of the pri- 
mary and of the secondary unstable modes. We could not 
achieve a reliable decomposition of the two m=2 modes be- 
cause of the noise and the nonlinear effects. Nevertheless, it is 
possible to estimate which of the two m =2 modes is dominant 
in a particular region of the disc. 

The unstable modes developing in A^-body simulations 
agree qualitatively with the linear analysis not only in the 
growth rate and the pattern speed, but also in their spatial ap- 
pearance. Figure |8] shows the contour plots of the surface den- 
sity for the 777=2, 3 and 4 Fourier components as determined 
from A^-body simulations (top frames) compared to the sur- 
face density contour plots for the unstable modes calculated by 
a linear stability analysis. For comparison, only the secondary 
777=2 linear mode is shown in the bottom left frame of Fig- 
ure [8^. As one can see, the spatial range, and winding of the 
unstable modes agree qualitatively in both A^-body and linear 
predictions. 

Figure [HJi shows a comparison of spatial distributions of 
perturbations within the central half -kiloparsec of the disc. The 
density distribution for the central bar-mode agrees with the 
analytically predicted density distribution plotted in the lower 
left panel. An agreement for density distributions of 777=3, and 
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Fig. 5. The dependence on time of m =2 (crosses), m =3 (triangles), and m =4 (squares) amplitudes of the global modes 
A, measured at different radii. In the central regions of the disk the fastest growing is the m =2 global bar-mode. At the disc 
periphery the three armed perturbation prevails over other competitors. The time evolution for the phases of m =2 (crosses) and 
m =3 (triangles) Fourier components 6 is shown in the lower frames. 



m=4 modes is also qualitatively acceptable but the details are 
different because of the limited resolution of the A^-body model 
discussed above. 



4. Summary 

We use high resolution A'^-body simulations to follow the dy- 
namics of growing spiral perturbations that develop in a colli- 
sionless disc from the initial noise perturbations. At least ten 




Fig. 8. a) Top frames - the surface density distribution for m=2, 3 and 4 Fourier components determined from A^-body simulations 
as compared to the surface density contour plots for the unstable modes determined from a linear stability analysis (bottom 
frames). For m=2 the secondary mode is shown (lower left), b) Density distributions within a central half-kiloparsec for the 
Fourier components shown in Fig.|8^. Here the primary m=2 mode is plotted at the lower left for comparison. 
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Fig. 6. it Upper and middle panels: growth rates and pattern 
speeds for m=2, m=3 and rn=4 Fourier components measured 
at different radii. In the central region of the disc, the growth 
rate and pattern speed of the m=2 Fourier component agree 
with the linear analysis prediction of the bar mode (soUd lines). 
Analytic predictions of the growth rates and pattern speeds are 
indicated at larger radii for m=2 (dashed lines), rn=3 (dotted 
Unes) and to=4 (dot-dashed lines) Fourier components. Lower 
panel: Points - the amplitudes of m=2, 3 and 4 - Fourier com- 
ponents as a function of radius taken from A'^-body simulations 
at time 40 Myr. The ampUtude profile for the most unstable 
m=2 central bar-mode agrees qualitatively well with the linear 
analysis (solid line). 
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Fig. 7. The growth rate and the pattern speed of the principal 
m =2 bar-mode as a function of a number of particles in N- 
body simulations. 



million particles with a minimum grid resolution of 128"'' are 
needed to reach a robust accuracy level so that the evolution of 
perturbations is no longer dominated by noise. 

Comparison of A'^-body simulations with the results of lin- 
ear stability analysis shows an agreement between both ap- 
proaches. The most unstable global bar-mode developing from 
the noise perturbations has a pattern speed and a growth rate as 
predicted by theory. Other unstable modes (the primary m=3, 
m=4 modes, and a weak secondary m=2 mode) quahtatively 
agree with the theoretical results. However the pattern speeds 
and growth rates of the secondary two-armed spiral, the three- 
armed and four-armed spirals are in lesser agreement with the 
theoretical values due to a stiU low resolution of the disc dy- 
namics in outer regions. 

We have demonstrated that a fast particle-mesh code Su- 
PERBOX is useful to follow a detailed dynamics of a collision- 
less stellar disc. Due to the large number of particles which 
can be used in SUPERBOX simulations, the noise level, and 
the numerical heating can be reduced to an insignificant level 
which allows us to foUow the complex dynamics of a coUision- 
less stellar disc. 

Theoretical models of JH are two dimensional, and have a 

finite mean rotation at the disc centre due to a strong positive 
gradient of azimuthal velocity dispersion. Stability properties 
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of the disc strongly depend on the inner boundary conditions of 
the disc. We plan to explore new theoretical models that have 
more realistic equilibrium properties. 
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